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Tlie parallel mutation-selection evolutionary dynamics, in which mutation and replication are 
independent events, is solved exactly in the case that the Malthusian fitnesses associated to the 
genomes are described by the Random Energy Model (REM) and by a ferromagnetic version of the 
REM. The solution method uses the mapping of the evolutionary dynamics into a quantum Ising 
chain in a transverse field and the Suzuki- Trotter formalism to calculate the transition probabilities 
between configurations at different times. We find that in the case of the REM landscape the 
dynamics can exhibit three distinct regimes: pure diffusion or stasis for short times, depending 
on the fitness of the initial configuration, and a spin-glass regime for large times. The dynamic 
transition between these dynamical regimes is marked by discontinuities in the mean-fitness as well 
as in the overlap with the initial reference sequence. The relaxation to equilibrium is described by 
an inverse time decay. In the ferromagnetic REM, we find in addition to these three regimes, a 
ferromagnetic regime where the overlap and the mean-fitness are frozen. In this case, the system 
relaxes to equilibrium in a finite time. The relevance of our results to information processing aspects 
of evolution is discussed. 

PACS numbers: 87.10.-e, 87.15. A-, 87.23. Kg, 02.50.-r 



I. INTRODUCTION 

Evolution on complex fitness landscapes has been ad- 
vanced as a key idea in recent mathematical approaches 
to evolution theory P, Concepts such as neutral 
networks and punctuated equilibrium, which are central 
to such disparate areas as molecular evolution and pa- 
leontology, can be brought together within that unify- 
ing framework. Since frustration and quenched disor- 
der combined together yield an almost infallible recipe 
to generate complexity [3| , the evolutionary dynamics on 
rugged fitness landscapes became a research topic in the 
statistical mechanics of disordered systems. 

There are a few reasons to consider spin-glass or ran- 
dom fitness landscapes in the study of molecular evolu- 
tion. Fitness functions are related to the binding affin- 
ity of a molecular replicator (a RNA-like molecule) to 
a non-specific replicase Q]- Given our present incapac- 
ity to predict affinity - single amino acid replacements 
may prevent binding altogether, increase affinity by or- 
ders of magnitude, or simply leave it unaffected [5'| ~ 
the assignment of fitness values chosen at random from 
some probability distribution seems to be the least biased 
course to introduce fitness in evolution models. In addi- 
tion, evolution in any fitness fitness landscape character- 
ized by a finite correlation length will resemble evolution 
in a random landscape when viewed at an appropriate 
coarse-grained scale of the sequence configuration space 
0. Finally, the analysis of the evolution on rugged or 
random fitness landscapes has produced dynamical pat- 
terns, such as punctuated equilibria (see [3,Ja]), that are 
actually observed in microbial populations Q . 

Building on a mapping between the infinite-population 



quasispecies model [10[ and an anisotropic two- 
dimensional Ising spin model in which the time t is one of 
the lattice dimensions 111, [13] , the evolutio nary version 
of Derrida's Random Energy Model (REM) [H, [li| was 
solved exactly in the infinite-time, equilibrium regime 
[isl . [l6| . Two distinct phases were found, corresponding 
to the selective and non-selective regimes that charac- 
terize a model that exhibits an error threshold transition 
(see also [l3|)- More recently, the dynamics of this model 
was investigated in the limit of strong selection using ap- 
proximate techniques [H, [H, [13, [Hj ■ Since strong selec- 
tion is not the only biologically relevant situation, and 
the selective regime is not the only important dynamic 
regime, a more general approach to the dynamics of the 
evolutionary version of REM is necessary. 

In this contribution, we explore a different mapping 
between evolutionary dynamics and statistical physics, 
namely, the mapping between the parallel mutation- 
selection scheme and the quantum Ising chain in a trans- 
verse field [2^, to solve exactly the dynamics on REM- 
like fitness landscape for all range of the selection and 
mutation parameters. This approach was already suc- 
cessfully used in the analysis of a simpler fitness land- 
scape, the Single-Peak fitness landscape [2^ and is based 
on the results of Refs. [13, [l^- A more recent applica- 
tion was the solution of the evolutionary dynamics in the 
case of symmetric fitness landscapes (i.e., the fitness is a 
function of the Hamming distance from a given reference 
sequence) ^26(1. 

Here we consider two landscapes, the ordinary REM 
landscape and the ferromagnetic REM landscape, where 
one of the energy levels of the ordinary REM is selected 
and changed to an energy value lower than the typical 
ground state energy. For both landscapes, we find that 
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the dynamic behavior for short times depends on the 
(Malthusian) fitness of the initial configuration. When 
the fitness of the initial configuration is higher than the 
mutation rate, the dynamics freezes at the initial con- 
figuration, resulting in a pattern of stasis (provided that 
the initial fitness is not the global maximum). Other- 
wise, when the initial fitness is lower than the mutation 
rate, the dynamics is characterized by a regime of pure 
diffusion in the sequence space. In the ordinary REM, 
we find that these short-time regimes - diffusion or sta- 
sis - change abruptly to a spin-glass-like regime which is 
associated to the equilibrium frozen phase of the REM. 
This dynamic transition is signaled by a discontinuity of 
the mean fitness as well as of the average overlap with 
the initial, reference sequence. Most importantly, we find 
that these quantities tend to their equilibrium values as 
1/t for large time t. In the case of the ferromagnetic 
REM we find up to four distinct dynamic regimes: dif- 
fusion or stasis, spin glass and ferromagnetic. As before, 
the transitions between these regimes are signaled by dis- 
continuities in the biologically relevant observables, but 
the system relaxes to the equilibrium ferromagnetic state 
in a finite time. 

Our results are interesting for the statistical physics as- 
pects of information theory as well [13, [H, [H H IM HI , 
as REM statistical physics gives a simple derivation of 
most of information theory results. The idea of coding 
via statistical mechanics is to construct a spin Hamilto- 
nian, which has a known ground state for a specific choice 
of the deterministic spin couplings (the non-trivial aspect 
of the problem is that the ground state of the Hamilto- 
nian should be robust to some degree of noise in those 
couplings). The ground state of the Hamiltonian could 
then be recovered from the starting configuration using 
some update dynamics, as in the case of associative mem- 
ory neural networks [33|. Thus by solving an evolution 
model we get as a by-product an analytical decoding dy- 
namics for optimal codes. 

The rest of the paper is organized as follows. In Sect. 
im we introduce the evolution equations for the parallel 
mutation-selection scheme and discuss its relation with 
the quantum Ising chain in a transverse field. The basic 
equations for the dynamics obtained using the Suzuki- 
Trotter formalism 34] arc introduced also in that sec- 
tion. In Sect, [m] we review the main results obtained in 
the analysis of the Single-Peak landscape f23l| as they un- 
derlie most of the arguments used in the solution of the 
more complex landscapes. In Sects . ITVl and IVl we present 
the exact solution of the evolutionary dynamics of the 
ordinary REM and of ferromagnetic REM, respectively. 
Finally, in Sect. IVII we summarize our main results and 
present some concluding remarks. 



II. ISING QUANTUM CHAIN FORMULATION 

In this contribution we consider the so-called parallel 
mutation-selection scheme in which mutation and selec- 



tion are considered as independent events [33, [3a, l37| . 
i.e., mutations can occur at any time during the exis- 
tence of a sequence, not only at the moment of replication 
as assumed in Eigen's molecular quasispecies model 1^ . 
As usual, we represent a molecule or genome of length 
by a sequence of binary digits (spins) sj, = ±1 with 
k — 1, . . . , A so that there are 2^ distinct molecules S"^ = 
{s\, . . . , sJy). In the parallel mutation-selection scheme, 
the relative frequencies of molecules i = 1, 
given by [s^ 



2^ are 



dt 



2" 



(1) 



where are Malthusian fitnesses, which can take on pos- 



itive as well as negative values 



and 



is the muta- 



tion rate from 5' to . Since mutations can connect only 
nearest neighboring sequences in the 2^-dimensional se- 
quence space, we choose = 7 if d[S^,S^) = 1; 
rriii = — A^7 and my = 0, otherwise. Here d[S'^,S^) 
is the Hamming distance between sequences S'' and 



and 7 is the mutation rate per site. As 



0, the 



- 1 for 
deter- 



dynamics ([IJ maintains the normalization X^jP; 
all t. Finally, we recall that = f (^s\, . . . , s" 
mines the so-called fitness landscape. 

A key observation at this stage is the finding that the 
non-linear dynamic system ([T]) can be reduced to a linear 
system 



where Hij = Hji = 



dxi \ ^ 



(2) 



using the transformation 



Xi {t) = pt (t) exp 



dr pj (r) 



(3) 



In practice, we solve the linear system ^ and then obtain 
the original sequence frequencies via the normalization 
Pi = Xij Xj. From the numerical perspective, solution 
of this linear system is straightforward, the sole limitation 
being the exponential increase in the number of equations 
with the sequence length A. 

For certain fitness landscapes, however, the case of in- 
finite length sequences can be solved analytically thanks 
to a cunning observation by Baake et al. who re- 

alized that the linear system Q can be mapped into an 
Ising quantum chain in a transverse magnetic field with 
spin interactions that depend on the specific choice of 
the fitness landscape. More pointedly, the linear system 
([5|) is equivalent to the evolution of the quantum system 
described by the Hamiltonian [22 1 



N 



^-7 E^^ 



A^ 



'N) 



(4) 
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where a^'^ stands for the Pauh spin operators acting in 
site k, i.e., ct^'"' = 1 (g) . . . 1 (g) a^'^^ (g) 1 . . . (g) 1 with ct^-^ in 
the A;th place. 

Introducing the time evolution operator T (t) — 
exp {—TLt) we can write a formal expression for the orig- 
inal molecules frequencies, namely, 



(5) 



be also the most studied fitness landscape in the quasis- 
pecies literature.) In this case, there is a single config- 
uration - the master sequence 5*" - with a high fitness 
value iVJo, whereas all other configurations have their 
fitness values set to zero. By choosing the master se- 
quence as S'^ = (f , . . . , I), the spin interactions for the 
SP landscape can be written as [2^ 



i=l 



fsp (cti, . . . , fJ 



NJo 



(9) 



{S^\Tit)\S^) andAA=E,,^ 



E,,^,, it)p^ (0) 



where Zji (t) 

guarantees the correct normahzation. Here \S^) = IxD® 
. . . (g Ixjv) where \x%) is an eigenstate of erf. 

Since there is an equivalence between the quantum 
Ising chain, such as that described by the Hamiltonian 
and a classical anisotropic two-dimensional Ising 
spin model [s?], then Baake et al.'s observation shows 
that there is a mapping between the parallel mutation- 
selection evolution scheme and the two-dimensional Ising 
model. It is interesting that a similar result holds for 
Eigen's quasispecies model as well [ll|, Ell • 

The challenge here is to calculate Zji (t) which is pro- 
portional to the probability that state 15') transitions 
to state 15^) in a time interval of length t. Henceforth 
we will refer to Zji as the transition amplitude between 
those states. In the case that the spins interaction, i.e., 
the term / (erf, . . . , cr^) in Eq. (jH), can be neglected we 
have T — > Tdi// = exp [7X^1 (^f ~ 1)] so that Zji can be 
readily evaluated [U (see also [23}), 

{S^ \%.,ff (t) \S^) = exp [N (m, t) - ^t)] (6) 



in the limit p ^ 00. What makes the SP problem ana- 
lytically solvable was the remarkable finding that Z can 
be written in a factorized form 

Z = V / dh{S^\%nt{t-h)\S'){S^Td,}f{ti)\S'). 

_ (10) 

We refer the reader to the Appendix of Ref. [23] for the 
detailed derivation of this result, which is based on the 
Trotter-Suzuki scheme .34;] (see also [13, [11]) that intro- 
duce infinitely many intermediate time steps between the 
initial and the final configurations. The key point for the 
factorization is the large p-spin interaction [see Eq. ([9])] 
in the Hamiltonian, which results in a ground-state con- 
figuration with fitness much greater than the fitness of 
typical configurations. This is obvious for the SP land- 
scape (except for the master, all configurations have zero 
fitness), but also holds for the REM landscape for which 
typical configurations have fitness on the order of iV^/^ 
whereas the ground-state has fitness on order of iV. 



where 

1 "t" 777, 1 — TTL 

(j) {m, t) = — ^ — In cosh {'^t) H ^ — In sinh {'yt) (7) 



and m is the overlap between configurations 5* and , 

he Ci 



i.e., m = J2k^k^k/-^- other hand, in the case 

the interactions are dominant we have T 
exp [—N-ft + f {af, . . . , aff) t] and so 



exp 



-N-ft + f {S') t] 4 



(8) 



For a given ensemble of initial configuration 5* our aim 
is to determine which configurations maximize the 
transition amplitude Zji (t). From Eq. ^ we can al- 
ready realize that the answer will depend only on the 
overlap between these two configurations and so the 
problem is reduced to finding the overlap m that maxi- 
mizes Zij. Henceforth we will refer to this maximum as 
Z = maxj Zji, thus omitting, for the sake of simplicity, 
the dependence on the initial configuration index i. Of 
course, because of the large N limit we have Z = Zji 
as well. 

The simplest fitness landscape for which the transition 
amplitudes Zji can be calculated exactly is the so-called 
Single-Peak (SP) fitness landscape. (The SP happens to 



III. SP FITNESS LANDSCAPE 

It is instructive to present the results for Z in the case 
of the single-peak landscape defined in the previous sec- 
tion and investigated in Ref. [23I ]. In particular, here 
we focus on the time dependence of the average overlap 
m between the configurations at time t and the initial 
configuration, a quantity which was not explored in that 
seminal work. First, we replace the sum over the final 
configurations by an integral over all possible values 
of the overlap between the final and the initial configu- 
rations, taking into account that for large N there are 
p (m) = exp [Nh (m)] distinct configurations for a fixed 
overlap m, where 

,, , I+m, 1-1-777 1 — 777, 1 — 777 

h(m) = — m — — In — - — . II 

^ ' 2 2 2 2 ^ ' 

Then we use Eqs. dH) and ^ to rewrite Z as 



Z = J dm J dtiexp[NFsp {m,ti 



)] (12) 



where 



Fsp = h{m) + (f> (777, h) --/t + Joit- h) 



(13) 



4 



The integrals over m and ti can be easily carried out 
for large N using Laplace's method as only the contribu- 
tion of the maximum of Fs p is relevant for the evaluation 
of Z. We find a maximum at the extreme of the ti inte- 
gration interval, i.e., ti — t, which, according to Eq. (|10p . 
corresponds to a regime of pure diffusion in the sequence 
space. In fact, maximization of Fgp with respect to m 
for ti =t yields 



m — exp (— 27t) , 



(14) 



from where we obtain Fsp = 0. Next, we consider the 
maximum within the integration intervals. The condition 
of maximum with respect to ti yields 



(1 + m) tanh {'^ti) + 



1 — TO 

tanh (7ti) 



-^^0. (15) 

7 



This is a quadratic equation for the unknown tanh (7ii) 
which has real solutions provided that Jq > j regard- 
less of the value of the other unknown, to. This is 
an evidence that this solution describes the selective 
regime of the parallel selection-mutation evolution model. 
Now, maximization of Fsp with respect to m yields 
TO = exp(— 27<i). The problem is that inserting this 
expression into Eq. ([TS]) results in Jq — 7, and so for this 
particular situation to can take on any value, whereas ti 
is given by Eq. ([T3]) and Fsp =0. 

To understand what is happening here we recall that 
the average overlap between the sequences in the qua- 
sispecies distribution at equilibrium and the master se- 
quence equals exactly 1 in the N ^ 00 limit. This is so 
despite the fact that the master sequences comprise only 
the fraction 1 — 7/J0 of the total number of sequences 
at equilibrium. In the same vein, the overlap to between 
the initial configuration and the equilibrium configura- 
tions become identical to the overlap TOq between the 
initial configuration and the master sequence, which is 
an arbitrary quantity defined by the initial conditions of 
the system. This explains the fact that the overlap to in 
Eq. (|15p is not specified by the maximization conditions 
for Jq = j: it is determined by the initial conditions, 

771 ~ TOq. 

To take into account the possibility that the dynamics 
reaches the close neighborhood of the master sequence, 
and so the overlap to (t) freezes at the value toq, in a fi- 
nite time we now calculate the transition amplitude Z^i 
between and the master sequence 5° in the time inter- 
val t, given that toq = J2k ^k^t/-^- Since the final state 
is fixed the entropic term (jlip must be dropped and we 
find 

Zo,(xe^p[N{^{mo,h)--ft + Jo{t-h))] (16) 

where ti is given by Eq. (|15p with to = toq. Recalling 
that Fsp — for the diffusive regime pi)) . the selective 
regime takes over at time t such that (j3(mo,ti) — + 
Joit-ti) > Figure [TJ which exhibits the time- 

dependence of the overlap to, summarizes these results 
for a particular choice of J0/7 and toq. The jump of the 




FIG. 1: Average overlap with the initial configuration as func- 
tion of the scaled time for J0/7 ~ 2, mo = 0.5 and (broken 
lines from bottom to top at jtdf ~ 0.763) N = 8, 12, 16 and 
24. The thick solid line is the theoretical prediction. 



overlap to at "/t^f signals the transition between the diffu- 
sive and the selective (frozen) dynamic regimes. We note 
that for the value of J0/7 = 2 used in the figure, the over- 
lap jumps up if Too > 0.112 and jumps down otherwise. 
The transition between the two regimes is continuous for 
Too ~ 0.122. As illustrated in Fig. [H the results of the 
numerical integration of the system of equations ([1]) using 
the Runge-Kutta method show a clear trend to converge 
to the theoretical predictions as the sequence length N 
increases. 

To conclude this brief overview of the parallel 
mutation-selection dynamics for the SP landscape we 
mention that the mean fitness R of the population is 
zero in the diffusive regime and i? = Jq — 7 (i.e., Jo 
times the frequency of master sequences in the popula- 
tion, 1 — 7/ Jo) in the selective phase. 

IV. REM FITNESS LANDSCAPE 



In this case the fitness landscape is given by [I3, uA 



/(si,. 



, sn) 



(17) 



where the couplings Ji-^,,,ip are Gaussian distributed ran- 
dom variables of zero mean and variance (jf^ ^ ) = 
J2p!/ (2iVP-i). Taking the limit p -> cx) in Eq. 
(fT7)) results in 2^ independent energy levels, E = 
— f (si, . . . , Sn), distributed by a Gaussian distribution 



w{E) = 



: exp 



(18) 



The equilibrium statistical mechanics of the quantum 
Ising model in a transverse field, Eq. with spin in- 
teractions given by Eq. ITTl) was studied in Ref. [2^ . 
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In the zero-temperature limit, which is the hmit rele- 
vant to our analysis, there are two phases, namely, a 
spin-glass or frozen phase that occurs for large J, and 
a paramagnetic phase. The discontinuous transition be- 
tween these phases takes place at J/7 — 1 / \/ln 2 « 1.201 
[2^ . Clearly, within our evolutionary interpretation of 
the model, this phase transition corresponds to the er- 
ror threshold phenomenon: for mutation rates greater 
than JVln2 the adaptive information stored in the fitness 
landscape no longer affects the frequencies of sequences 
in the population, which become uniformly distributed. 

Because of the presence of the quenched random vari- 
ables Jii...ij,, the derivation of the equations for the dy- 
namics of the REM is more involved compared with that 
for the SP landscape. A rigorous derivation of the factor- 
ization (1101) can be done using the Trotter-Suzuki scheme 
[H, m, [SJI as in the SP fitness landscape [2^ . Here we 
give a qualitative derivation based on the general results 
described in Sect. [TTl 

Effectively, we simply must replace / (S**) in Eq. ([5]) 
by Ei (i.e., by the energy of the initial configuration S"*), 
which amounts to replacing Jq in Eq. by Ei/N. 

Since our final results must be averaged over the ener- 
gies of the final configurations [see Eq. (ITU)) ] we have 



where 




h (to) 



E 



X exp 



E^ 
' N.P 



E{t-h) 



X exp [Nh (m) + Ncj) (to, ti) - Njt] . 
Here the theta function enforces the constraint 



h (to) 



(- 



> 



(19) 



(20) 



which guarantees that the average number of configura- 
tions with energy E and overlap to with the initial con- 
figuration, given by exp (^Nh{m) — E^/NJ^), is exponen- 
tially large, so that Z becomes a self-averaging quantity 
[1^ . Hence Eq. (|19p describes a typical situation 
of the dynamics at a fixed time t. Strictly, this equa- 
tion is valid when the energy of the initial configuration 
Ei is larger then —7. We will discuss later the correc- 
tions necessary to describe the case where this condition 
is violated. In addition, we assume that the overlap toq 
between the ground-state configuration and the initial 
configuration is zero. 

The crucial step now is to evaluate the integral over 
E in Eq. ([TO]) via Laplace's integration for fixed to and 
ti. Noting that the result of the integration depends on 
whether the critical point E* = NJ^ {t — ti) /2 satisfies 
or not the constraint (|20[) and dropping trivial multiplica- 
tive factors, we rewrite Z as 

^D = J dti J dm exp[NFD {m,ti)] (21) 



FD = h (to) + (to, h) --ft + J^{t- hf /4 (22) 



provided that 



/i(to)- j2(i-ii)2/4>0, 



(23) 



ZsG^ I dh f dm exp[NFsG{m,h)] (24) 
Jo J -1 



where 



FsG^(bim,ti)-jt + Jy/h^{t-h) (25) 
in the case that 



h (to) 



hf/A < 0. 



(26) 



Equation (|24p results when the critical point E* is out- 
side the energy integration interval, and so the argu- 
ment of the exponential is maximized by the energy E at 
the extreme of that interval, namely, E — NJ^Jh (to). 
As before, for large N we need to find the values of to 
and ti that maximize the arguments of the exponentials, 
Eqs. ([m and (|25p . The correct solution is then the 
one that corresponds to the largest of Z^) and Zsq, i.e., 
Z = max{Z£), Zsg]- 

We begin with the analysis of Zu-, Eq. (|2ip . Calcula- 
tion of the extreme of Fd (to, t\) with respect to to and 
tx yields to = exp {-2^t + 47^/72) and t\ = t - 
so that Fo — — 7^/J^ < 0. Next we must evaluate Fo 
at the upper extreme of the ti integration interval, i.e., 
at ti = t. This is the pure diffusion regime discussed 
in Sect. Iml which results in Eq. ^ and Fd = 0. The 
lower extreme = yields F^ — > — 00 and the extremes 
of the TO integration interval (i.e., m — ±1) need not 
be considered because /i (to = ±1) = and so the con- 
dition (|23p is violated. We note that the solution given 
by Eq. p4p , which describes the pure drift or diffusion in 
the sequence space, exists for all parameter values since 
h (to) > and so the inequality (P5)) is always satisfied. 
In addition, since solution (fTl]) yields the largest value of 
the exponent Fd, the other solutions must be discarded. 

We turn now to the analysis of Zsg, Eq. ([M)) . As 
before, we start by the maximization of Fsg (™j ^i) with 
respect to both integration variables, to and ti. At the 
maximum, we find that the values of these variables are 
given by the solution of the equations 

Intanh {^t,) + i In ( ^-±^\ = (27) 



2 



and 



(1 + to) tanh (7ti) -|- 



1 — TO 



1 — TO 



2J 



tanh(7ii) 7 



yjh (m) = 0. (28) 



These equations have to be solved numerically, but the 
solution is simple because Eq. (pS)) can be rewritten as 
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a quadratic equation y = tanh (7^1) < 1 which then can 
be written expHcitly in terms of the unknown m. Re- 
gardless of the value of to, this quadrat ic eq uation has 
real solutions provided that J/7 > l/-\/ln2 and so we 
identify this dynamic regime with the frozen spin-glass 
phase of the quantum version of REM [24| . In addition, 
in the case that both roots of y are physical (i.e., less 
than 1), our numerical analysis indicates that we should 
always choose the smaller root since it corresponds to the 
largest value of the exponent Fsg- 

To conclude the analysis of Zsg, we must consider the 
contributions from the extremes of the integration inter- 
vals. The extreme ti = t is discarded because it violates 
condition (p6)) . whereas the contribution of ii = can 
be ignored because it yields Fsg ~^ ~oo. Regarding 
the extremes to = ±1, we find that in this case -Fsg is 
maximum when ti takes on its extreme value, ti — t. 
This corresponds to the contribution from the border 
h (to) — J^(t— ii)^/4 = which we will discuss in detail in 
the Appendix [X] Our numerical analysis indicates, how- 
ever, that the border contribution can be ignored since 
it yields an exponent Fb [see Eq. (IA2[) ] which is always 
smaller than the exponents obtained using the solution 
of the Eqs. ^ and (^5)) . 

In addition to the average overlap to between the initial 
and the configuration at time t, we can calculate the 
time dependence of the mean fitness R of the sequence 
population as well. The reasoning to derive R is sketched 
as follows. For the diffusive regime we have R — since 
the average fitness of any large sample of configurations 
visited in this regime is clearly zero for the REM fitness 
landscape. To estimate the mean fitness in the selective 
regime we just note the equivalence between the results 
for the single-peak landscape [see Eqs. ([T^ and 
and for the selective phase [see Eqs. p4)) and ((25)) ] if we 
identify Jq with an effective, time-dependent single-peak 
fitness value Jeff — J\/h (m). (Note that for i ^ 00 we 
have TO — > so that Je/ / —^ Jv'ln2, which is the ground- 
state fitness value of the REM.) Since the population is 
formed by master copies with fitness value NJ^f / as well 
as by clouds of mutants with much smaller fitness (on 
the order of N^^'^) the mean fitness of the population in 
the selective regime becomes 

R = Jy/h{m) - J, (29) 

in accord with the well-known result for the (parallel) 
version of the single-peak landscape. 

In summary, for fixed J/7 and "ft we must solve the 
saddle-point equations (|27|1 and ([28]) to obtain the expo- 
nent Fsg (as well the saddle-point equation (|A3P given 
in Appendix but we have already mentioned that its 
contribution must be discarded) and then compare with 
the exponent of the diffusive regime Fjj — 0. If Fsg > 
we pick the value of to associated to the selective regime, 
otherwise we pick the diffusion solution given by Eq. p4p . 



E 




FIG. 2: Average overlap with the initial configuration as func- 
tion of the scaled time "/t for (thin solid lines from top to 
bottom at "/t = 2) J/7 = 2, 3, 4 and 5. The thick solid line 
is the function exp {—2"/t) which describes the overlap in the 
diffusive regime, J/7 < 1/ Vln 2. 



A. Analysis of the results 

Figure[2]illustrates the typical time evolution of the av- 
erage overlap m with the initial configuration. We recall 
that the fitness of this initial configuration must be less 
than 7 and its overlap with the ground-state configura- 
tion must be zero. As expected, for small "ft the diffusive 
regime dominates and so to is given by Eq. (|14p . As "ft 
increases further, the selective regime takes over rather 
abruptly, as shown by the discontinuity of the overlap 
TO at a critical time value "ftds- This bizarre behavior, 
which occurs also in the single-peak fitness landscape, is 
consequence of our characterization of the dynamics in a 
very large-dimensional sequence space by a single param- 
eter: no such discontinuous behavior is observed when 
following the time evolution of the individual sequence 
frequencies, pi for i = 1, . . . , 2^ . 

Since the properly scaled critical time "ftds at which the 
discontinuity of the overlap m (and, consequently, of the 
mean fitness R) takes place can be used to separate the 
regions of validity of the two distinct dynamic regimes, 
in Fig. [3] we present the dynamic 'phase diagram' of the 
parallel evolutionary version of REM. As expected, "ftds 
diverges as the J/ 7 approaches the value 1 / Vln 2 which 
yields the equilibrium phase boundary between the para- 
magnetic and the frozen spin glass phases. For large J/7 
we find that "ftds vanishes as {J/"f)~ ■ Also of interest is 
the size of the overlap jump at "ftds, shown in Fig. [J] The 
fact that this quantity exhibits a maximum could already 
be inferred from Fig. [51 since the overlap to tends to 1 
or in both dynamic regimes when J/7 approaches its 
extreme values. 
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FIG. 3: Scaled critical time 'ytds at which the discontinuous 
dynamical transition between the diffusive and the selective 
regimes takes place as function of the dimensionless parameter 
J/7. For J/7 < l/\/ln 2 only the diffusive regime occurs. The 
selective regime is dominant in the region t > tds (i.e., above 
the solid line). 



FIG. 5: Scaled mean fitness R/'y of the REM landscape as 
function of the scaled time for J/7 = 4 and (broken lines 
from top to bottom at -ytds « 0.292) iV = 12, 16, 20 and 24. 
The thick solid line is the theoretical prediction. 
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FIG. 4: Size of the overlap discontinuity Am at t = tds as 
function of the dimensionless parameter J/7. The maximum 
of this curve occurs at J/7 « 2.504. For for large J/7 we find 
that Am vanishes as (J/'j)"^. 



B. Numerical integration 

To complement our theoretical analysis, which is exact 
for infinite sequence lengths, we have carried the direct 
numerical integration of the linear system of ordinary 
equations ^ for sequence lengths up to iV = 24 using 
the fourth-order Runge-Kutta integrator [3§|. The sta- 
bility of the numerical procedure benefited greatly from 
the fact the differential equations are linear. For < 10 
we can find all eigenvectors and eigenvalues of the sym- 
metric matrix H and so solve the dynamics exactly for 



Yt 

FIG. 6: Average overlap with the initial configuration as func- 
tion of the scaled time 'yt for J/7 = 4 and (broken lines from 
bottom to top at jtds ~ 0.292) iV = 12, 16, 20 and 24. The 
thick solid line is the theoretical prediction. 



any t within an arbitrarily high numerical precision. Of 
course, the two numerical methods yield identical results 
provided that is not too large. Figures [5] and [6] sum- 
marize our numerical results for the J/7 = 4. For each 
time t the data in these figures represent the average 
over 10* independent samples. The samples differ by the 
fitness values assigned to each configuration. For all sam- 
ples the initial configuration was such as to have fitness 
value less than 7 and zero overlap with the ground-state 
configuration. 

Figure [5] is reassuring because the crossings of the lines 
for distinct N indicate the onset of a threshold phe- 
nomenon in the thermodynamic limit N 00. In par- 
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FIG. 7: Values of 'yt at which the mean fitnesses of sequences 
of length A'' and A'' + 2 intersect shown as function of 
for N = 12, 14, 16, 18, 20 and 22 for J/7 = 4. The symbols A 
and O identify the first and the second crossings, respectively 
(see Fig.[SJ, whereas the filled symbols indicate the theoretical 
predictions. The solid lines are the linear fittings used to 
obtain the extrapolated values at 1/N = 0. 



ticular, to reproduce the analytical predictions, the first 
intersection point should tend to 7t = whereas the 
second should tend to ^tds ~ 0.292. To verify whether 
our data exhibit the correct trend, we show in Fig. [7] 
the values of at which the mean fitness curves inter- 
sect for successive values of TV. The extrapolation to 
N 00 yields •yt — — O.OliO.Ol for the first crossing and 
ytds = 0.30 ± 0.01 for the second crossing. The agree- 
ment with the theoretical prediction is excellent, given 
the short sequence lengths used in the numerical integra- 
tion. Oddly enough, the dependence of the overlap m 
on the sequence length N, shown in Fig. [SI does not ex- 
hibit the characteristic crossings that signalize the onset 
of a threshold phenomenon in the thermodynamic limit, 
although the curve for TV = 24 already begins to take a 
shape that resembles the theoretical prediction. It seems 
that much larger sequence lengths are needed in order 
we can obtain clear evidence of a threshold phenomenon 
using the overlap data. We refer the reader to Ref. [4l| 
for a full analysis of the finite size effects of the error 
threshold transition of the quasispecies model. 



C. High-fitness initial configuration 

To complete our analysis of the calculation of Z we 
consider now the situation in which the energy of the 
initial configuration Ei is such that Ei < —7, i.e., this 
configuration has a relatively high fitness. In this case, 
we find a new paramagnetic (but non-diffusive) regime 
where the system stays in the original configuration with 



a probability proportional to [see Eq. ([8])] 

Z = exp [N i~E, ~ 7) t] (30) 

and mean fitness R — —Ei. Since the argument of the ex- 
ponential is always positive, this new regime replaces the 
diffusive regime altogether. As t increases, it eventually 
becomes replaced by the selective regime at some thresh- 
old time t'ds > tds- However, the dependence oit'^^ on the 
particular value Ei makes this case rather unattractive, 
as compared with the case where the initial configuration 
has a low-fitness value. 

Finally, we note that the probability that Ei < —7 is 

ierfc ^7/ J\/]vj which tends to 1/2 for large N, hence we 

could more simply distinguish the two situations - high 
and low initial fitness - by verifying whether the energy 
of the initial configuration is positive or negative. More 
importantly, these two cases are equally likely and our 
penchant for the low-fitness initial configuration here is 
justified only by the generality of the results obtained in 
that case. 



D. Relaxation to equilibrium 

The approach to the equilibrium state as 7i — > cxd is 
particularly interesting because it is related to the speed 
of evolution, i.e., how long it takes for a random se- 
quence to reach the global maximum of a rugged fitness 
landscape. In contrast with the finite-time dynamics de- 
scribed before, the results of the asymptotic analysis do 
not depend on the specific value of the fitness of the initial 
configuration. We still require, however, that the initial 
configuration has zero overlap with the ground state. 

As we focus on the limit of large jt, the relevant equa- 
tions to describe the system dynamics are Eqs. (P7|) and 
(PS)) . From Fig. [5] we can see that m — > in this limit 
and so Eq. yields 

y = K — \/ — \ (31) 

where y — tanh (7ii) and k — J^/h\2/^. Since y < 1 for 
K > 1, ^1 is finite and then, taking the limit < ^ 00 in 
Eq. dUl), we find 

/ ln21ny\ 1 

™^(, — —)Tf ^''^ 

This important result indicates that relaxation to equilib- 
rium, which is characterized by the ground-state config- 
uration together with a cloud of very close mutant con- 
figurations (the quasispecies distribution) is given by a 
power law with exponent —1. 

An estimate of the speed of evolution can be obtained 
by considering the prefactor of 1/7^ in Eq. ((32|) . We 
find that this prefactor vanishes at the extremes k = 1 
and K ^ 00, and reaches a maximum at k w 1.810 or 
J/7 w 2.174. (If we had included the curve for, say. 
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FIG. 8: Asymptotic dependence of the average overlap on 'yt 
for J/7 = 4 and (broken lines from bottom to top) A'^ = 6, 8 
and 10. The thick solid line is the theoretical prediction, 
which yields m oc 1/ {-yt) [see Eq. . 



J/7 = 1.5 in Fig. [5] we could have observed this non- 
monotonic behavior there.) This parameter setting cor- 
responds then to the slowest convergence to equilibrium, 
i.e. the minimum speed of evolution. The maximum 
speed is obtained by setting k — > cx) (or J/7 — > 00) which 
amounts to taking a vanishingly small mutation rate. 

To check whether a similar power-law scaling for large 
times holds also for (infinite) populations of finite length 
sequences, we solved the linear system ([2]) through the 
direct diagonalization of H, which is feasible only for 
relatively small sequence lengths. Figure [8] summarizes 
our numerical results, which represent the average over 
10'' independent samples. We find that the finite N 
data is very well fitted by the scaling law m oc (7t)~"" 
with ag = 1.77, as — 1-36 and aio = 1.17. As- 
suming that aoo = 1 we find that these exponents, 
in turn, are described perfectly by the function ajq = 
1 -I- 7.48 exp (— 0.38A^), which indicates a very rapid ap- 
proach to the value of the infinite-length exponent. 



V. FERROMAGNETIC REM LANDSCAPE 

In the ferromagnetic REM [l^, [23| we choose a partic- 
ular configuration, say S*" = (1, . . . , 1), and set its fitness 
value to JqN . The other 2^ — 1 configurations are as- 
signed random fitness values —E with E distributed by 
the Gaussian distribution (jlSp . From the evolutionary 
modeling perspective, the new fitness level produces a 
gap in the fitness landscape, which, as we show here, re- 
sults in nontrivial dynamic consequences. 

The quantum spin version of the ferromagnetic REM is 



defined by the Hamiltonian ([4]) with the fitness function 
/ (si, . . . , Sat) = Jii...ipSi-i ■ . ■ Si^ 

il <Z2---<'ip 



(33) 



where the multispin couplings Ji-^.,,ip are defined as in 
Sect. IIVI This fitness landscape is thus a linear combi- 
nation of the SP and REM landscapes. The equilibrium 
statistical mechanics of the quantum ferromagnetic REM 
was studied in Ref. [25], where the condition for the ex- 
istence of the ferromagnetic phase at zero-temperature 
(i.e., for S'^ be the ground-state configuration) was found 
to be Jo > JVln 2. Here we will consider only parameter 
settings that satisfy this condition. 

As in the previous cases, we use the decomposition of 
the transition amplitude Z, Eq. (jlOp . to solve the dynam- 
ics for the overlap m between the initial and the typical 
configurations at time t. The important change is that 
now the sum over the final configurations S-' in Eq. (jlOp 
does not include the master sequence S'^ , which must be 
considered separately. Hence we find that Z is given by a 
sum of two terms, the first is the REM contribution, Eq. 
P^ . and the second is the SP contribution, Eq. (fTBl) . In 
particular, we focus on the case toq = only, so that the 
latter equation becomes 



Zp = exp 



^ l^^sinh(27ti) 
2 2 



jt + Jo{t-t[) 



(34) 

where y' — tanh {-ft'-^ ) is given by the quadratic equation 
{y'f - 2 Jo2/77 + 1 = (see Eq. ^ with m = mo = 0), 
which has real solutions for Jq > 7. As the interesting 
situation is one where the spin-glass solution, give n by 
Eqs. ^ and exists as weU, so that J/7 > l/\/hi2, 
we have 



^>:^Vh^>i, 

7 7 



(35) 



so the existence of the ferromagnetic and spin-glass 
phases guarantees that y' is real. 

To obtain the time evolution of m (t) we first calculate 
In ZsG using Eq. ([24]) and In Zp using Eq. ([34]) for fixed t. 
(We take logarithms here because the relevant quantities 
are the expressions in the arguments of the exponentials 
that define the transition amplitudes.) If these two quan- 
tities arc negative, wc choose the diffusive solution, Eq. 
P^ . If \nZp > In ZgQ then the system is in the ferro- 
magnetic regime and so 771 = 0; otherwise we choose m 
given by the spin-glass solution, Eqs. (P7)l and 

As the setting Jo > JVhi 2 implies that the equi- 
librium phase is the ferromagnetic one, one must have 
TO = for large 'yt. On the other hand, if the fitness 
of the initial configuration is not greater than 7 (which 
we tacitly assume in this section), the diffusive regime 
dominates for small 'yt. The question is then whether 
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FIG. 9: Average overlap m as function of scaled time jt. The 
solid curve for Jo /7 — J/7 = 3 shows a transition between the 
diffusive and the spin-glass regimes at yt = 0.538, and a tran- 
sition between the spin-glass and the ferromagnetic regimes 
at 'yt = 0.820. The broken curve for J0/7 = 3 and J/j — 2 
shows a situation where there is a direct transition between 
the diffusive and ferromagnetic regimes at 'yt = 0.693. 
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FIG. 10: The spin-glass regime interfaces the diffusive and the 
ferromagnetic regimes only in the region of parameters located 
below the thick solid curve. This curve begins at the point 
(l/\/ln2, 1) and diverges as J^ for large J. The thin solid 
straight line is Jo — JVln 2, below which the ferromagnetic 
phase is absent. 



an intermediate, spin-glass regime appears between these 
extremes. The answer is given in Fig. |9l which indicates 
that the appearance of the intermediate regime depends 
on the values of the parameters J and Jq. To determine 
the region in the space of parameters (J/7, Jq/j) where 
the spin-glass regime interfaces the other two regimes, we 
have to calculate the value of J0/7 such that the time t^s 
at which the transition from the diffusive to the spin-glass 
regime coincides with the time t^f at which the diffusive 
regime transitions to the ferromagnetic one. Note that 
tds depends only on J as shown in Fig. 01 As for the 
time t = tdf at which the transition between the diffusive 
and the ferromagnetic regimes takes place, it can be cal- 
culated analytically by setting \nZp — (see [11]). The 
final result is 



1 



4 {ko - 1) 



kq In 



1 



In [4 (4 1)] 



(36) 



where kq = Jo/l- The procedure for searching the values 
of J0/7, for fixed J/7, at which t^s = tdf is implemented 
numerically and the result is shown in Fig. 1101 Above 
the thick solid line there are only two dynamic regimes, 
the diffusive and the ferromagnetic, and the time tdf at 
which the transition occurs is given by Eq. ()36|) . In what 
follows we will concentrate on the study of the dynamics 
for the parameters in the region below that line, where 
the three dynamic regimes are present. In particular we 
will focus on the transition between the spin-glass and the 
ferromagnetic regimes, which happens at time t — tgf. 

Before we offer an analytical approximation to tgf it is 
instructive to study numerically its dependence on J0/7 
and J/7. This is shown in Fig. 1101 from where it be- 
comes clear that tsf is defined in a narrow region of the 




FIG. 11: Scaled time at which the spin-glass transitions to 
the ferromagnetic for (left to right) J/7 = 2,3 and 4. This 
transition occurs only within a limited region of the parameter 
space, as illustrated in Fig. 1101 The divergences occur at 
Jo = J\/In2. 



parameter space, determined by the conditions that the 
ferromagnetic phase exists and that the spin-glass regime 
interfaces the diffusive and ferromagnetic regimes. Fig- 
ure[TT]shows the discontinuity of the overlap at tsf - Since 
the overlap in the ferromagnetic regime is zero, the jump 
Am is actually the overlap in the spin-glass phase. 

The divergence of t^f and the vanishing of Am = m 
as Jo approaches JVln 2 allows us to derive an analytical 
expression for t^f in this limit. In fact, Eq. p2p already 
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Of course, since e — > we must set ^ 1 in the former 
equation. Keeping leading order terms in 6 = kq — 1, Eq. 
([55)1 becomes ^tdf ^ (In 2) /S, whereas Eq. (^0]) reduces 

1/2 

to jtsf ^ {S hi 2/e) ' . Equating these resuhs yields S — 

1/3 

(e In 2) so that the minimum time to reach the optimal 
fitness situation is 

(In 2)^/^ , , 

7imm = ^ ^1/3 ■ (41) 

This expression is valid only in the limits J0/7 — > 1 and 
J/7 ^ 1/Vhr2. 

VI. CONCLUSION 



FIG. 12: The overlap discontinuity Am at t = tsf for (left 
to right) J/7 = 2, 3 and 4. Since the overlap is zero in the 
ferromagnetic regime, Am equal the overlap in the spin-glass 
regime. 



provides an explicit expression for m, namely, 

Vhi21ny 



Jt 



(37) 



since is/ is large. The equation that defines t^f is ob- 
tained by equating InZp/N to Fsg, 



In 



sinh {2-ft[) 



sinh (27ti) 



-f m In tanh {■yti) + 2J ^ h {rn)ti — 
+2Jot[ - 2 (jo - J^h (m)) tsf. (38) 



Recalling that for e = (^ Jq — Jvhi2j / Jq — > we have 
K ^ kq and so ^ ii , we rewrite this expression as 



t 



In tanh {'yti) 
' 2Jn^ 



(39) 



where we used h (m) ^ In 2 for m — )■ 0. Finally, inserting 
m from Eq. (|37p into this expression yields 



2e 



■ In Ko 



(40) 



As m ~ 1/ {'^tgf) we find the typical mean-field result 
m ~ e^/^ at the transition. 

Since tsf (or tc;/, depending on the parameter settings) 
is the waiting time for evolution to lead the system close 
to its optimal fitness situation, it is interesting to see 
whether this waiting time can be minimized by a proper 
choice of parameters (the mutation rate, for example). 
For fixed kq , we note that tsf must satisfy the constraint 
tsf > tdf in the case the spin-glass regime is present. 
Hence the minimum waiting time tmm is obtained by 
equating the waiting times given in Eqs. (|36p and (HO]). 



Most of the techniques from statistical mechanics em- 
ployed in the study of evolutionary models, such as 
Eigen's quasispecies model, are manageable only in the 
stationary regime t — > 00 (see, e.g., [H, [H, [2^ lioj). 
Although the equilibrium analysis provides valuable in- 
sights into the behavior of these models, a complete study 
of the dynamics is indispensable as evolution is all about 
species dynamics, after all. 

In this contribution we present an exact solution for 
the evolutionary dynamics in an extremely rugged fit- 
ness landscape, Derrida's Random Energy Model (REM) 
[1^. The evolutionary model studied is the quasis- 
pecies model with a parallel mutation-selection scheme, 
in which mutations are decoupled from replication [36^. 
This scheme can be mapped in Ising quantum chain in 
a transverse field [l^], and the dynamics can be solved 
using the Suzuki- Trotter formalism as done in the case 
of the Single-Peak (SP) landscape [1^. In fact, the sim- 
ilarity between the SP and REM-like fitness landscapes 
regarding their steady-state distributions - they are iden- 
tical within the accuracy ~ 1/ \/N - is well-known [l^ , 
and here we explore it to derive the evolutionary dynam- 
ics on the REM landscape using the SP landscape results. 

The (infinite) population is initially homogeneous, i.e., 
all sequences are identical to a reference sequence chosen 
such that its overlap with the highest-fitness sequence is 
zero. In addition, most of our results are based on the 
assumption that the fitness of this reference sequence is 
negative. We note that in the parallel mutation-selection 
scheme, we have a Malthusian fitness which basically 
measures the difference between the reproduction and 
death rates, and so can take on positive and negative 
values as well. 

At each time t, the population is characterized by the 
average overlap with the reference sequence m (t) as well 
as by the mean fitness R (t) . As expected, in the case the 
initial configuration has low fitness (i.e., the fitness value 
is less than the mutation rate per site 7), the dynamics 
for small t corresponds to a random drift in the sequence 
space with m decreasing exponentially with increasing t 
[see Eq. (HI])]. Selection, which encodes information in 
the fitness landscape, has no role in the diffusive regime. 
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We find, quite remarkably, that m undergoes a discon- 
tinuous transition at some finite t = tj^g (see Fig. [2]) when 
the dynamics enters a spin-glass regime in which m van- 
ishes as 1/t for large t. As opposed to the SP landscape 
(and to the ferromagnetic version of REM; see below), 
the dynamics needs an infinite time to reach the regions 
close to optimal fitness sequence. When the initial con- 
figuration already has high fitness (i.e., the fitness value 
is greater than 7) the diffusive regime is replaced by a 
pattern of stasis: the dynamics freezes at the initial con- 
figuration (i.e., ni — 1 and R — —Ei) for a certain length 
of time ti > tds where ti — ti {Ei) and then undergoes a 
discontinuous transition to the spin-glass regime. 

In addition to the REM fitness landscape, we con- 
sidered also the somewhat more realistic ferromagnetic- 
REM landscape which, as it is clear from Eq. can 
be viewed as a simple combination of the REM and SP 
fitness landscapes. For some parameter settings (see Fig. 
|9]), we find three distinct dynamic regimes: the diffusive, 
spin-glass and the ferromagnetic regimes. The transitions 
between these regimes are signaled by discontinuities of 
the overlap as well as of the mean fitness. In a parame- 
ter setting such that the equilibrium phase is the ferro- 
magnetic one, the time to reach the optimal sequence is 
finite, as in the SP case, but diverges near the (equilib- 
rium) transition points. As in the case of the ordinary 
REM, the diffusive regime is replaced by stasis when the 
initial configuration has a high fitness value. 

The discontinuous transitions between the different dy- 
namic regimes are similar to the punctuations, during 
which evolution proceeds very rapidly, observed in finite 
population simulations 0,3- In fact, one of the first the- 
oretical models to reproduce the punctuated equilibrium 
phenomenon made explicit use of the effect of random 
genetic drift, which results from the finitude of the pop- 
ulation, to pro mote the transition between alternative fit- 
ness peaks [i^l . Since our results were derived within the 
infinite-population assumption, this process cannot be 
responsible for the observed punctuations. Alternatively, 
punctuations are predicted by models in which initially 
low frequency beneficial mutation becomes dominant in 
a few generations after a certain frequency threshold is 
overcome [4^ 13] . This is the process responsible for the 
punctuations observed in our model as well as in micro- 
bial population experiments [§] . 

The infinite size population assumption behind the 
quasispecies-like evolution model considered here is a the- 
oretical approximation only, and finite population size 
effects are undoubtedly important. The discrete-time 
evolutionary dynamics on a REM-like fitness landscape 
has been extensively investigated in the literature for 
the finite population case [3, [l^j [13, [U, EE]- Ana- 
lytical approximations and numerical simulations have 
yielded many interesting results about record statistics 
and crossover transitions. We note, however, that the 
exact solution of the deterministic model exhibits a much 
richer dynamical structure. It would be interesting to see 
whether there are any vestiges of the discontinuous tran- 



sitions in the case of finite but large populations and, in 
particular, how the coalescent time statistics are affected 
by these regime changes j45| . 

While our main interest is the investigation of the evo- 
lutionary dynamics, our results bear on information the- 
ory as well d^, as they can be viewed as the exact 
analytical solution for the decoding process (relaxation 
to the ferroma gnet ic configuration) of optimal codes 
(27I m, m, I30I l3ll|. In particular, we conjecture that 
Eq. (|¥T|) is universal for some classes of dynamics near 
the error threshold-like transitions. In fact, the connec- 
tion of evolution models with information theory was first 
pointed out by Eigen, who actually used information the- 
oretical arguments to derive an expression for the error 
threshold in the SP landscapes [I^]- More recently, the 
relation between molecular biology and information the- 
ory was discussed in Refs. [47ll48|. 

In general, almost any fitness landscape can be qual- 
itatively identified with one of three classes: ferromag- 
netic, spin-glass and ferromagnetic spin-glass like. The 
first class, which exhibits a finite relaxation time to the 
optimum, seems too simplistic to bear on real biologi- 
cal situations. The spin-glass fitness landscape, on the 
other hand, exhibits an infinite relaxation time to the 
optimum, which then could never be reached by the evo- 
lutionary dynamics. The third class, which combines the 
complexity of the spin-glass landscape with a finite re- 
laxation time, seems to be the preferable one from the 
evolutionary perspective. Thus natural selection seems 
to choose the type of fitness landscape that works more 
efficiently as an information processing system. 

The picture that emerges from computer experiments 
with digital organisms [43| resembles the case of SG fit- 
ness case. Although this random macroevolution sce- 
nario may be described by a spin-glass fitness landscape, 
Nature's preference seems to be for the ferromagnetic 
spin-glass landscape, as manifested, for example, by pro- 
tein evolution. In fact, it is known that proteins dif- 
fer substantially from random heteropolymers, and that 
random heteropolymers can be described by the ordi- 
nary REM, whereas biological polymers are described 
quite well by the ferromagnetic REM [s^, [5l|. Hence 
the genome, which codes the information to assemble 
the proteins, reveals ferromagnetic or ferromagnetic spin- 
glass like fitness landscape. This is close to the idea of 
channels in evolution [52l |. During the evolution, there 
are large rearrangements of the genome, in addition to 
the point substitutions considered here. This large trans- 
positions resemble the multi-scale optimization in com- 
puter science: perhaps Nature takes advantage of large 
gene rearrangements, whenever the search for new opti- 
mum by means of simple substitutions becomes too slow 
[5^. These large events, as well as the simultaneous point 
mutations in two or three adjacent sites, are permitted 
because they practically do not affect the error threshold, 
while the evolution dynamics changes drastically, e.g., a 
relaxation time of about 10^ years in the case of sin- 
gle point mutations is reduced to 100 years when triplet 
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adjacent mutations are allowed [53|. We note that si- 
multaneous mutations in two or three random sites yield 
the same slow relaxation as in the case of a single point 
substitution. 

In the theory of computation, optimization problems 
are classified as polynomially solvable if the relaxation 
time to the optimum scales with some power of the prob- 
lem size, and as NP-complete or as NP-hard otherwise, 
i.e., when the relaxation time scales exponentially with 
the problem size [13] • Hence the channel-like evolution 
schemes are not only ferromagnetic-type fitness, but also 
resemble the fast computational schemes of polynomial 
class. 
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APPENDIX A: BORDER CONTRIBUTION TO Z 

To complete the analysis of the integrals in Eqs. ((2T|) 
and p4|) for large N, we must consider the contribution 
from the border h (to) = {t — ti) /4. Clearly, in this 
case we have Zjy = Zsg = Zb with 



Zb = J dm exp [NFb (r 



where 



2h{m) + 0(TO,ti) - "ft 



(Al) 



(A2) 



and ti — ti (to) is a function of m given by the border 
equation. Maximization of Fb with respect to m yields 



7ln ft™ 



In 



to 



— lntanh(7ii) = 



(1 + to) tanh (7ti) 



tanh (7ii) 



(A3) 



which must be solved numerically. At the extremes 
TO = ±1 we have ti = t and so Fb = In [1 ± exp {—2^t)] — 
In 2 < 0. Our extensive numerical analysis comparing 
the three exponents Fd = 0, Fsg and Fb indicates that 
whenever Fb > we have Fsg > Fb and so the contri- 
bution from the border can be neglected in comparison 
with those from the inner saddle points discussed in the 
main text. 



